1 Import library

Code

# data viz
library(ggplot2)
library(ggh4x)
library(ggOceanMaps)
library(patchwork)
library(viridis)
library(ggpubr)
library(grid)
library(viridis)
library(ggnewscale)

# table
library(gt)

# map
library(rnaturalearth)
library(ggmap)
library(ggsn)
library(sf)
library(sp)
library(smoothr)

# kernel density
library(eks)
library(ggdensity)

# stat
library(Hmisc)
library(circular)
library(CircStats)

# char
library(stringr)

# solar angle
library(oce)

# data wrangling
library(tidyr)
library(dplyr)
library(purrr)
library(forcats)
library(lubridate)

# model
library(lme4)
library(fitdistrplus)

2 Setting up custom function

2.1 `windrose`

Code

windrose <-
  function(data_to_plot,
           grid = NULL,
           set_title = NULL,
           legend_position = "none",
           facet = F) {
    # this code comes from Roxanne
    uniqhours <- 1:24 * (360 / 24)

    # trick to align hours om the graph
    data_to_plot <- rbind(data_to_plot[2:nrow(data_to_plot), ], data_to_plot[1, ])

    for (i in 1:24) {
      # turn hours to radians
      if (i == 1) {
        temp <- rep(uniqhours[i], data_to_plot$nb_ind_hour[i])
        day_night <- rep("night", data_to_plot$nb_ind_hour[i])
      } else {
        temp <- c(temp, rep(uniqhours[i], data_to_plot$nb_ind_hour[i]))
        day_night <- c(day_night, rep(
          if_else(between(i, 7, 20), "day", "night"),
          data_to_plot$nb_ind_hour[i]
        ))
      }
    }
    data2 <- data.frame(direction = temp)

    deg <- 15 # choose bin size (degrees/bin)
    dir.breaks <- seq(0 - (deg / 2), 360 + (deg / 2), deg) # define the range of each bin

    # assign each direction to a bin range
    dir.binned <-
      cut(data2$direction,
        breaks = dir.breaks,
        ordered_result = TRUE
      )
    # generate pretty lables
    dir.labels <- as.character(c(seq(0, 360 - deg, by = deg), 0))

    # replace ranges with pretty bin lables
    levels(dir.binned) <- dir.labels

    # Assign bin names to the original data set
    data2$dir.binned <- dir.binned

    # add origin if any
    if (facet) {
      data2$origin <- unique(data_to_plot$origin)
    }

    # set up max value
    maxvalue <- 35

    # initialise the plot
    plt.dirrose_2 <- ggplot()

    # check if grid
    if (!is.null(grid)) {
      plt.dirrose_2 <- plt.dirrose_2 +
        geom_hline(
          yintercept = grid,
          colour = "grey20",
          linewidth = .2
        )
    }
    plt.dirrose_2 <- plt.dirrose_2 +
      geom_vline(
        xintercept = c(seq(1, 24, 2)),
        colour = "grey30",
        linewidth = 0.2
      ) + # 24 vertical lines at center of the 30? ranges.
      geom_hline(
        yintercept = maxvalue,
        colour = "white",
        linewidth = .5
      ) + # Darker horizontal line as the top border (max).
      # On top of everything we place the histogram bars.
      geom_bar(
        data = data2,
        aes(x = dir.binned, fill = day_night),
        width = 1,
        colour = "black",
        linewidth = 0.3
      ) +
      # geom_bar(data = data2, aes(x = dir.binned2), width = 1, colour="black", size = 0.3,fill="salmon",alpha=0.9) +
      scale_x_discrete(
        drop = FALSE,
        labels = c(
          0, "", 2, "", 4, "", 6, "", 8, "", 10, "", 12, "", 14, "", 16, "", 18, "", 20, "", 22, ""
        )
      ) +
      scale_fill_manual(values = c("white", "darkgrey"), name = "Time of day") +
      labs(x = "Time (hours)", y = "Count", title = set_title) +
      coord_polar(start = -(deg / 2) * (pi / 180))

    # if facet
    if (facet) {
      plt.dirrose_2 <- plt.dirrose_2 +
        facet_wrap2(. ~ origin)
    }

    # Wraps the histogram into a windrose
    plt.dirrose_2 <- plt.dirrose_2 +
      theme_bw() +
      theme(
        legend.position = legend_position,
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        panel.background = element_blank(),
        legend.key.size = unit(0.5, "cm")
      )

    # return
    return(plt.dirrose_2)
  }

2.2 `getNightLength`

Based on here: https://stackoverflow.com/questions/59566553/calculate-night-length-in-r

Code

getNightLength <- function(survey_date, latitude, longitude) {
  # survey_date <- strptime(survey_date, "%d/%m/%y")
  time_origin <- strptime("01/01/2000 12:00", "%d/%m/%Y %H:%M")
  n <- as.numeric(difftime(survey_date, time_origin, units = "days"))
  Jstar <- n - longitude / 360
  M <- (357.5291 + 0.98560028 * Jstar) %% 360
  C <- 1.9148 * sin(M * 2 * pi / 360) +
    0.0200 * sin(2 * M * 2 * pi / 360) +
    0.0003 * sin(3 * M * 2 * pi / 360)
  lambda <- (M + C + 180 + 102.9372) %% 360
  Jtransit <- as.double(2451545.000) +
    as.double(Jstar) +
    as.double(0.0053 * sin(M * 2 * pi / 360)) -
    as.double(0.0069 * sin(lambda * 4 * pi / 360))
  sindelta <- sin(lambda * 2 * pi / 360) * sin(23.44 * 2 * pi / 360)
  delta <- asin(sindelta) * 360 / (2 * pi)
  cos_omega <- (sin(-0.83 * 2 * pi / 360) - (sindelta * sin(latitude * 2 * pi / 360))) /
    (cos(latitude * 2 * pi / 360) * cos(delta * 2 * pi / 360))
  omega <- acos(cos_omega) * 360 / (2 * pi)
  Jrise <- Jtransit - omega / 360
  Jset <- Jtransit + omega / 360
  return((1 - (Jset - Jrise)) * 24)
}

3 Import data

Let’s load data_dive, i.e. the output of data_wrangling.qmd, and filter only on animals leaving from Ano Nuevo.

Code

# import the data
data_dive <- readRDS("../export/data_dive.rds")

# filter on seals departing from ANNU
data_dive <- data_dive %>%
  filter(DepartureLocation == "ANNU")

Code

# get solar angle when data location
solar_angle_inter <- data_dive %>%
  filter(!is.na(Lat)) %>%
  mutate(solar_angle = sunAngle(date, Long, Lat)$altitude) %>%
  select(id, DiveNumber, solar_angle)

# merge that with data_dive
data_dive <- merge(data_dive,
  solar_angle_inter,
  by = c("id", "DiveNumber"),
  all.x = T
)

# add day-night
# https://sciencepickle.com/earth-systems/sun-earth-connection/declination-circles/sunrise-sunset-and-twilight/
data_dive <- data_dive %>%
  mutate(day_night = if_else(solar_angle < -18, "night", "day"))

Let’s add a condition based on the difference between the maximum depth reached and the bathymetry to define a benthic dive, i.e. if this difference is above xxx m then is not considered as a benthic dive.

Code

data_dive <- data_dive %>%
  # change DiveType to include condition on difference
  # between max depth and bathymetry (<50m)
  mutate(
    DiveType_50 = data.table::fifelse(DiveType == 3 &
      (abs(bathy) - Maxdepth) > 50, 0, DiveType),
    DiveTypeName_50 = data.table::fifelse(
      DiveTypeName == "Benthic" &
        (abs(bathy) - Maxdepth) > 50,
      "Transit",
      DiveTypeName
    ),
    # same with <100m
    DiveType_100 = data.table::fifelse(DiveType == 3 &
      (abs(bathy) - Maxdepth) > 100, 0, DiveType),
    DiveTypeName_100 = data.table::fifelse(
      DiveTypeName == "Benthic" &
        (abs(bathy) - Maxdepth) > 100,
      "Transit",
      DiveTypeName
    ),
    # same with <150m
    DiveType_150 = data.table::fifelse(DiveType == 3 &
      (abs(bathy) - Maxdepth) > 150, 0, DiveType),
    DiveTypeName_150 = data.table::fifelse(
      DiveTypeName == "Benthic" &
        (abs(bathy) - Maxdepth) > 150,
      "Transit",
      DiveTypeName
    ),
    # same with <200m
    DiveType_200 = data.table::fifelse(DiveType == 3 &
      (abs(bathy) - Maxdepth) > 200, 0, DiveType),
    DiveTypeName_200 = data.table::fifelse(
      DiveTypeName == "Benthic" &
        (abs(bathy) - Maxdepth) > 200,
      "Transit",
      DiveTypeName
    )
  )

4 Maps

Code

# check if background_ggoceanmaps exist
if (file.exists("../export/background_ggoceanmap.rds")) {
  trip <- readRDS("../export/background_ggoceanmap.rds")
} else {
  # using ggOceanMaps
  trip <- basemap(
    limits = c(170, -110, 20, 59),
    bathymetry = TRUE,
    shapefiles = "Arctic",
    rotate = TRUE,
    grid.col = NA
  )

  # Make the graticules:
  lims <- attributes(trip)$limits
  graticule <- sf::st_graticule(
    c(lims[1], lims[3], lims[2], lims[4]),
    crs = attributes(trip)$proj,
    lon = seq(-180, 180, 45),
    lat = seq(-90, 90, 10)
  )

  # Plot
  trip <- trip +
    geom_sf(data = graticule, color = "grey50")

  trip$layers <- trip$layers[c(1, 3, 2)]

  # save result
  saveRDS(trip, "../export/background_ggoceanmap.rds")
}

Code

# data use to compute kernel density estimation
df_kernel_dens <- data_dive %>%
  # only with location data
  filter(!is.na(Lat) & DiveTypeName == "Benthic") %>%
  # select only the required columns
  select(lat = Lat, long = Long, id, DiveTypeName)

Code

# transform data into sf object
df_kernel_dens_sf <- st_as_sf(
  df_kernel_dens,
  coords = c("long", "lat"),
  crs = st_crs(4326)
)

# kernel density estimation
df_kernel_dens_sf_kde <- st_kde(
  x = df_kernel_dens_sf,
  # boundary = ne_coastline(scale = 110, returnclass = "sf"),
  H = diag(c(
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 1]
    ),
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 2]
    )
  ) / 4)^2
)

# https://github.com/r-spatial/sf/issues/1762
sf::sf_use_s2(FALSE)

# plot
test <- trip +
  # new scale
  new_scale_fill() +
  # kernel
  # ggplot(data=ne_coastline(scale = 110, returnclass = "sf"))+geom_sf()+
  geom_sf(
    data = st_get_contour(
      # geospatial kernel
      df_kernel_dens_sf_kde,
      # kde_test,
      # probabilities
      cont = c(50, 80, 95, 99)
    ),
    # display
    aes(fill = label_percent(contlabel)),
    alpha = 0.7
  ) +
  # same colour bar
  scale_fill_viridis_d(option = "plasma") +
  # legend
  labs(fill = "Probs") +
  # no display alpha
  guides(alpha = "none") +
  # title
  labs(
    title = "All Benthic Dives",
    subtitle = paste(nrow(df_kernel_dens), "dives")
  )

# reorder layers
test$layers <- test$layers[c(1, 2, 4, 3)]

# plot
test

Figure 1: Kernel density plots of the all benthic dives performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park.

Code

# data use to compute kernel density estimation
df_kernel_dens <- data_dive %>%
  # only with location data
  filter(!is.na(Lat) & DiveTypeName_50 == "Benthic") %>%
  # select only the required columns
  select(lat = Lat, long = Long, id, DiveTypeName)

Code

# transform data into sf object
df_kernel_dens_sf <- st_as_sf(
  df_kernel_dens,
  coords = c("long", "lat"),
  crs = st_crs(4326)
)

# kernel density estimation
df_kernel_dens_sf_kde <- st_kde(df_kernel_dens_sf,
  H = diag(c(
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 1]
    ),
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 2]
    )
  ) / 4)^2
)

# https://github.com/r-spatial/sf/issues/1762
sf::sf_use_s2(FALSE)

# plot
test <- trip +
  # new scale
  new_scale_fill() +
  # kernel
  geom_sf(
    data = st_get_contour(
      # geospatial kernel
      df_kernel_dens_sf_kde,
      # probabilities
      cont = c(50, 80, 95, 99)
    ),
    # display
    aes(fill = label_percent(contlabel)),
    alpha = 0.7
  ) +
  # same colour bar
  scale_fill_viridis_d(option = "plasma") +
  # legend
  labs(fill = "Probs") +
  # no display alpha
  guides(alpha = "none") +
  # title
  labs(
    title = "Benthic Dives with bathy-max_depth<50",
    subtitle = paste(nrow(df_kernel_dens), "dives")
  )

# reorder layers
test$layers <- test$layers[c(1, 2, 4, 3)]

# plot
test

Figure 2: Kernel density plots of benthic dives with a difference between the bathymetry and the maximum depth of 50 m performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park.

Code

# data use to compute kernel density estimation
df_kernel_dens <- data_dive %>%
  # only with location data
  filter(!is.na(Lat) & DiveTypeName_100 == "Benthic") %>%
  # select only the required columns
  select(lat = Lat, long = Long, id, DiveTypeName)

Code

# transform data into sf object
df_kernel_dens_sf <- st_as_sf(
  df_kernel_dens,
  coords = c("long", "lat"),
  crs = st_crs(4326)
)

# kernel density estimation
df_kernel_dens_sf_kde <- st_kde(df_kernel_dens_sf,
  H = diag(c(
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 1]
    ),
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 2]
    )
  ) / 4)^2
)

# https://github.com/r-spatial/sf/issues/1762
sf::sf_use_s2(FALSE)

# plot
test <- trip +
  # new scale
  new_scale_fill() +
  # kernel
  geom_sf(
    data = st_get_contour(
      # geospatial kernel
      df_kernel_dens_sf_kde,
      # probabilities
      cont = c(50, 80, 95, 99)
    ),
    # display
    aes(fill = label_percent(contlabel)),
    alpha = 0.7
  ) +
  # same colour bar
  scale_fill_viridis_d(option = "plasma") +
  # legend
  labs(fill = "Probs") +
  # no display alpha
  guides(alpha = "none") +
  # title
  labs(
    title = "Benthic Dives with bathy-max_depth<100",
    subtitle = paste(nrow(df_kernel_dens), "dives")
  )

# reorder layers
test$layers <- test$layers[c(1, 2, 4, 3)]

# plot
test

Figure 3: Kernel density plots of benthic dives with a difference between the bathymetry and the maximum depth of 100 m performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park.

Code

# data use to compute kernel density estimation
df_kernel_dens <- data_dive %>%
  # only with location data
  filter(!is.na(Lat) & DiveTypeName_150 == "Benthic") %>%
  # select only the required columns
  select(lat = Lat, long = Long, id, DiveTypeName)

Code

# transform data into sf object
df_kernel_dens_sf <- st_as_sf(
  df_kernel_dens,
  coords = c("long", "lat"),
  crs = st_crs(4326)
)

# kernel density estimation
df_kernel_dens_sf_kde <- st_kde(df_kernel_dens_sf,
  H = diag(c(
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 1]
    ),
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 2]
    )
  ) / 4)^2
)

# https://github.com/r-spatial/sf/issues/1762
sf::sf_use_s2(FALSE)

# plot
test <- trip +
  # new scale
  new_scale_fill() +
  # kernel
  geom_sf(
    data = st_get_contour(
      # geospatial kernel
      df_kernel_dens_sf_kde,
      # probabilities
      cont = c(50, 80, 95, 99)
    ),
    # display
    aes(fill = label_percent(contlabel)),
    alpha = 0.7
  ) +
  # same colour bar
  scale_fill_viridis_d(option = "plasma") +
  # legend
  labs(fill = "Probs") +
  # no display alpha
  guides(alpha = "none") +
  # title
  labs(
    title = "Benthic Dives with bathy-max_depth<150",
    subtitle = paste(nrow(df_kernel_dens), "dives")
  )

# reorder layers
test$layers <- test$layers[c(1, 2, 4, 3)]

# plot
test

Figure 4: Kernel density plots of benthic dives with a difference between the bathymetry and the maximum depth of 150 m performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park.

Code

# data use to compute kernel density estimation
df_kernel_dens <- data_dive %>%
  # only with location data
  filter(!is.na(Lat) & DiveTypeName_200 == "Benthic") %>%
  # select only the required columns
  select(lat = Lat, long = Long, id, DiveTypeName)

Code

# transform data into sf object
df_kernel_dens_sf <- st_as_sf(
  df_kernel_dens,
  coords = c("long", "lat"),
  crs = st_crs(4326)
)

# kernel density estimation
df_kernel_dens_sf_kde <- st_kde(df_kernel_dens_sf,
  H = diag(c(
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 1]
    ),
    MASS::bandwidth.nrd(
      sf::st_coordinates(df_kernel_dens_sf)[, 2]
    )
  ) / 4)^2
)

# https://github.com/r-spatial/sf/issues/1762
sf::sf_use_s2(FALSE)

# plot
test <- trip +
  # new scale
  new_scale_fill() +
  # kernel
  geom_sf(
    data = st_get_contour(
      # geospatial kernel
      df_kernel_dens_sf_kde,
      # probabilities
      cont = c(50, 80, 95, 99)
    ),
    # display
    aes(fill = label_percent(contlabel)),
    alpha = 0.7
  ) +
  # same colour bar
  scale_fill_viridis_d(option = "plasma") +
  # legend
  labs(fill = "Probs") +
  # no display alpha
  guides(alpha = "none") +
  # title
  labs(
    title = "Benthic Dives with bathy-max_depth<200",
    subtitle = paste(nrow(df_kernel_dens), "dives")
  )


# reorder layers
test$layers <- test$layers[c(1, 2, 4, 3)]

# plot
test

Figure 5: Kernel density plots of benthic dives with a difference between the bathymetry and the maximum depth of 200 m performed by female northern elephant seals during their foraging trips separated by each of the four dive types. Level indicates the concentration of dives with red showing the highest concentration. Seals all departed and returned from Año Nuevo State Park.

5 Windrose

Code

# let's add a column with the local_time
data_windrose <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName == "Benthic")

# figure
windrose_benthic <- data_windrose %>%
  group_by(id) %>%
  filter(DiveTypeName == "Benthic") %>%
  mutate(time = as.numeric(format(date_tz, format = "%H"))) %>%
  group_by(time) %>%
  summarise(nb_ind_hour = n()) %>%
  mutate(origin = "Benthic") %>%
  windrose(., grid = c(3000, 6000, 9000), legend_position = "top", facet = T)

# display
windrose_benthic

Figure 6: Circular histogram plots displaying the time when they perform benthic dives across their foraging trip.

Code

# perform test
test_benthic <- r.test(conversion.circular(
  circular(
    data_windrose %>%
      group_by(id) %>%
      filter(DiveTypeName == "Benthic") %>%
      mutate(time = hour(date_tz) +
        minute(date_tz) / 60 +
        second(date_tz) / (60 * 60)) %>%
      pull(time),
    units = "hours"
  ),
  units = "radians"
))

# convert back into human reading hour
conv_test <- conversion.circular(
  circular(test_benthic$r.bar,
    units = "rad"
  ),
  units = "hours"
)
# print
print(paste0("The average time is ", round(conv_test[[1]], 1), " (p-value: ", round(test_benthic$p.value, 1), ")"))

[1] "The average time is 0.5 (p-value: 0)"

Code

# let's add a column with the local_time
data_windrose <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_50 == "Benthic")

# figure
windrose_benthic <- data_windrose %>%
  group_by(id) %>%
  filter(DiveTypeName == "Benthic") %>%
  mutate(time = as.numeric(format(date_tz, format = "%H"))) %>%
  group_by(time) %>%
  summarise(nb_ind_hour = n()) %>%
  mutate(origin = "Benthic") %>%
  windrose(., grid = c(1000, 2000, 3000), legend_position = "top", facet = T)

# display
windrose_benthic

Figure 7: Circular histogram plots displaying the time when they perform benthic dives (diff<50) across their foraging trip.

Code

# perform test
test_benthic <- r.test(conversion.circular(
  circular(
    data_windrose %>%
      group_by(id) %>%
      filter(DiveTypeName_50 == "Benthic") %>%
      mutate(time = hour(date_tz) +
        minute(date_tz) / 60 +
        second(date_tz) / (60 * 60)) %>%
      pull(time),
    units = "hours"
  ),
  units = "radians"
))

# convert back into human reading hour
conv_test <- conversion.circular(
  circular(test_benthic$r.bar,
    units = "rad"
  ),
  units = "hours"
)
# print
print(paste0("The average time is ", round(conv_test[[1]], 1), " (p-value: ", round(test_benthic$p.value, 1), ")"))

[1] "The average time is 0.1 (p-value: 0)"

Code

# let's add a column with the local_time
data_windrose <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_100 == "Benthic")

# figure
windrose_benthic <- data_windrose %>%
  group_by(id) %>%
  filter(DiveTypeName == "Benthic") %>%
  mutate(time = as.numeric(format(date_tz, format = "%H"))) %>%
  group_by(time) %>%
  summarise(nb_ind_hour = n()) %>%
  mutate(origin = "Benthic") %>%
  windrose(., grid = c(1000, 2000, 3000), legend_position = "top", facet = T)

# display
windrose_benthic

Figure 8: Circular histogram plots displaying the time when they perform benthic dives (diff<100) across their foraging trip.

Code

# perform test
test_benthic <- r.test(conversion.circular(
  circular(
    data_windrose %>%
      group_by(id) %>%
      filter(DiveTypeName_100 == "Benthic") %>%
      mutate(time = hour(date_tz) +
        minute(date_tz) / 60 +
        second(date_tz) / (60 * 60)) %>%
      pull(time),
    units = "hours"
  ),
  units = "radians"
))

# convert back into human reading hour
conv_test <- conversion.circular(
  circular(test_benthic$r.bar,
    units = "rad"
  ),
  units = "hours"
)
# print
print(paste0("The average time is ", round(conv_test[[1]], 1), " (p-value: ", round(test_benthic$p.value, 1), ")"))

[1] "The average time is 0.1 (p-value: 0)"

Code

# let's add a column with the local_time
data_windrose <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_150 == "Benthic")

# figure
windrose_benthic <- data_windrose %>%
  group_by(id) %>%
  filter(DiveTypeName == "Benthic") %>%
  mutate(time = as.numeric(format(date_tz, format = "%H"))) %>%
  group_by(time) %>%
  summarise(nb_ind_hour = n()) %>%
  mutate(origin = "Benthic") %>%
  windrose(., grid = c(1000, 2000, 3000), legend_position = "top", facet = T)

# display
windrose_benthic

Figure 9: Circular histogram plots displaying the time when they perform benthic dives (diff<150) across their foraging trip.

Code

# perform test
test_benthic <- r.test(conversion.circular(
  circular(
    data_windrose %>%
      group_by(id) %>%
      filter(DiveTypeName_150 == "Benthic") %>%
      mutate(time = hour(date_tz) +
        minute(date_tz) / 60 +
        second(date_tz) / (60 * 60)) %>%
      pull(time),
    units = "hours"
  ),
  units = "radians"
))

# convert back into human reading hour
conv_test <- conversion.circular(
  circular(test_benthic$r.bar,
    units = "rad"
  ),
  units = "hours"
)
# print
print(paste0("The average time is ", round(conv_test[[1]], 1), " (p-value: ", round(test_benthic$p.value, 1), ")"))

[1] "The average time is 0.1 (p-value: 0)"

Code

# let's add a column with the local_time
data_windrose <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_200 == "Benthic")

# figure
windrose_benthic <- data_windrose %>%
  group_by(id) %>%
  filter(DiveTypeName == "Benthic") %>%
  mutate(time = as.numeric(format(date_tz, format = "%H"))) %>%
  group_by(time) %>%
  summarise(nb_ind_hour = n()) %>%
  mutate(origin = "Benthic") %>%
  windrose(., grid = c(1000, 2000, 3000), legend_position = "top", facet = T)

# display
windrose_benthic

Figure 10: Circular histogram plots displaying the time when they perform benthic dives (diff<200) across their foraging trip.

Code

# perform test
test_benthic <- r.test(conversion.circular(
  circular(
    data_windrose %>%
      group_by(id) %>%
      filter(DiveTypeName_200 == "Benthic") %>%
      mutate(time = hour(date_tz) +
        minute(date_tz) / 60 +
        second(date_tz) / (60 * 60)) %>%
      pull(time),
    units = "hours"
  ),
  units = "radians"
))

# convert back into human reading hour
conv_test <- conversion.circular(
  circular(test_benthic$r.bar,
    units = "rad"
  ),
  units = "hours"
)
# print
print(paste0("The average time is ", round(conv_test[[1]], 1), " (p-value: ", round(test_benthic$p.value, 1), ")"))

[1] "The average time is 0.1 (p-value: 0)"

6 Boxplot of the benthic dive rate (day vs. night)

Since we thought of accoutering for the variation in dive duration to explore how benthic dive rate varies across day-time vs night-time, I just looked at the difference in dive duration between both period of the day.

Code

data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  ggplot(aes(x = day_night, y = Dduration)) +
  # the boxplot
  geom_boxplot(outlier.shape = NA) +
  # zoom
  coord_cartesian(ylim = c(0, 2600)) +
  # labels
  labs(x = "Time of the day", y = "Dive duration (s)")

Code

broom::tidy(wilcox.test(
  data_dive %>%
    # get rid of data without location information
    filter(!is.na(Lat)) %>%
    # at night
    filter(day_night == "night") %>%
    # get dive duration
    pull(Dduration),
  data_dive %>%
    # get rid of data without location information
    filter(!is.na(Lat)) %>%
    # at night
    filter(day_night == "day") %>%
    # get dive duration
    pull(Dduration)
))

Graph

Code

data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName == "Benthic") %>%
  mutate(
    hour = lubridate::hour(date),
    day = as.Date(date)
  ) %>%
  group_by(id, day, hour, day_night) %>%
  summarise(
    nb_benthic_dive = n(),
    mean_dive_duration = mean(Dduration),
    .groups = "drop"
  ) %>%
  group_by(id, day_night) %>%
  summarise(
    nb_average_benthic_dive = mean(nb_benthic_dive),
    mean_dive_duration = mean(mean_dive_duration),
    .groups = "drop"
  )


data_plot %>%
  ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +
  geom_boxplot(outlier.shape = NA) +
  labs(
    x = "Time of the day",
    y = "Average number of benthic dives per seal per hour"
  )

Figure 11: Boxplot of the number of (all) benthic dives per hour and individual accross day-time and night-time

Model

Code

# determine what distribution would best describe our data
descdist(data_plot$nb_average_benthic_dive, discrete = FALSE)

summary statistics
------
min:  1   max:  4.333333 
median:  1.555556 
mean:  1.669479 
estimated sd:  0.4221527 
estimated skewness:  1.473928 
estimated kurtosis:  6.414337

Figure 12: Result from `descdist` function to determine the best distribution to use for data modelling

Based on Figure 12, it seems that a gamma distribution would be the best to model this dataset.

\[ HourlyBenthicDiveRate \sim DayNight + MeanDiveDuration + (1|ID) \tag{1}\]

Based on a table shape like the following:

ID	HourlyBenthicDiveRate	DayNight	MeanDiveDuration
A	10	Night	20
A	5	Day	21
B	11	Night	20
B	4	Day	22

Code

# model
mdl_all <- glmer(
  nb_average_benthic_dive ~ day_night +
    scale(mean_dive_duration) +
    (1 | id),
  family = Gamma,
  data = data_plot %>% mutate(day_night = factor(day_night,
    # ordered = T,
    levels = c(
      "night",
      "day"
    )
  ))
)
# summary
broom.mixed::tidy(mdl_all) %>%
  gt() %>%
  tab_header(title = "GLMM")

effect	group	term	estimate	std.error	statistic	p.value
GLMM
fixed	NA	(Intercept)	0.64660385	0.005864671	110.254069	0.000000e+00
fixed	NA	day_nightday	-0.02527306	0.004559170	-5.543347	2.967447e-08
fixed	NA	scale(mean_dive_duration)	0.11434250	0.004412862	25.911191	4.981557e-148
ran_pars	id	sd__(Intercept)	0.06307474	NA	NA	NA
ran_pars	Residual	sd__Observation	0.11146813	NA	NA	NA

Go Further…

Code

data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName == "Benthic") %>%
  # I'm not so sure about this part..
  mutate(period_duration = if_else(
    day_night == "night",
    getNightLength(date, Lat, Long),
    24 - getNightLength(date, Lat, Long)
  )) %>%
  group_by(id, day_night) %>%
  summarise(
    nb_benthic_dive = n(),
    mean_dive_duration = mean(Dduration),
    mean_period_duration = mean(period_duration),
    .groups = "drop"
  )

\[ NbBenthicDive \sim DayNight + MeanDiveDuration + MeanDayNightDuration + (1|ID) \]

Based on a table shape like the following:

ID	NbBenthic	DayNight	MeanDiveDuration	MeanDayNightDuration
A	10	Night	20	7
A	5	Day	21	17
B	11	Night	20	6
B	4	Day	22	18

For more information: https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#overdispersion

Code

# model
mdl_all_extra <- glmer(
  nb_benthic_dive ~ day_night +
    scale(mean_dive_duration) +
    scale(mean_period_duration) +
    (1 | id),
  family = "poisson",
  data = data_plot %>% mutate(day_night = factor(day_night,
    # ordered = T,
    levels = c(
      "day",
      "night"
    )
  ))
)
# summary
broom.mixed::tidy(mdl_all_extra) %>%
  gt() %>%
  tab_header(title = "GLMM")

effect	group	term	estimate	std.error	statistic	p.value
GLMM
fixed	NA	(Intercept)	5.41274751	0.041968504	128.97166	0.000000e+00
fixed	NA	day_nightnight	-0.72679883	0.006652905	-109.24533	0.000000e+00
fixed	NA	scale(mean_dive_duration)	0.27370713	0.007918629	34.56497	8.494370e-262
fixed	NA	scale(mean_period_duration)	0.08335783	0.002524341	33.02162	3.975539e-239
ran_pars	id	sd__(Intercept)	0.83810984	NA	NA	NA

Code

qqnorm(residuals(mdl_all_extra))
qqline(residuals(mdl_all_extra))

Graph

Code

data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_50 == "Benthic") %>%
  mutate(
    hour = lubridate::hour(date),
    day = as.Date(date)
  ) %>%
  group_by(id, day, hour, day_night) %>%
  summarise(
    nb_benthic_dive = n(),
    mean_dive_duration = mean(Dduration),
    .groups = "drop"
  ) %>%
  group_by(id, day_night) %>%
  summarise(
    nb_average_benthic_dive = mean(nb_benthic_dive),
    mean_dive_duration = mean(mean_dive_duration),
    .groups = "drop"
  )


data_plot %>%
  ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +
  geom_boxplot(outlier.shape = NA) +
  labs(
    x = "Time of the day",
    y = "Average number of benthic dives per seal per hour"
  )

Figure 13: Boxplot of the number of (diff<50) benthic dives per hour and individual accross day-time and night-time

Model

Code

# determine what distribution would best describe our data
descdist(data_plot$nb_average_benthic_dive, discrete = FALSE)

summary statistics
------
min:  1   max:  5.166667 
median:  2.555556 
mean:  2.581992 
estimated sd:  0.6479778 
estimated skewness:  0.3352311 
estimated kurtosis:  3.477161

Figure 14: Result from `descdist` function to determine the best distribution to use for data modelling

Based on Figure 14, it seems that a gamma distribution would be the best to model this dataset.

\[ HourlyBenthicDiveRate \sim DayNight + MeanDiveDuration + (1|ID) \tag{2}\]

Based on a table shape like the following:

ID	HourlyBenthicDiveRate	DayNight	MeanDiveDuration
A	10	Night	20
A	5	Day	21
B	11	Night	20
B	4	Day	22

Code

# model
mdl_50 <- glmer(
  nb_average_benthic_dive ~ day_night +
    scale(mean_dive_duration) +
    (1 | id),
  family = Gamma,
  data = data_plot %>% mutate(day_night = factor(day_night,
    levels = c(
      "night",
      "day"
    )
  ))
)
# summary
broom.mixed::tidy(mdl_50) %>%
  gt() %>%
  tab_header(title = "GLMM")

effect	group	term	estimate	std.error	statistic	p.value
GLMM
fixed	NA	(Intercept)	0.41544372	0.003997844	103.916932	0.000000e+00
fixed	NA	day_nightday	-0.01594999	0.003819268	-4.176192	2.964297e-05
fixed	NA	scale(mean_dive_duration)	0.08109884	0.002975447	27.256018	1.410142e-163
ran_pars	id	sd__(Intercept)	0.03423312	NA	NA	NA
ran_pars	Residual	sd__Observation	0.14226941	NA	NA	NA

Go Further…

Code

data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_50 == "Benthic") %>%
  # I'm not so sure about this part..
  mutate(period_duration = if_else(
    day_night == "night",
    getNightLength(date, Lat, Long),
    24 - getNightLength(date, Lat, Long)
  )) %>%
  group_by(id, day_night) %>%
  summarise(
    nb_benthic_dive = n(),
    mean_dive_duration = mean(Dduration),
    mean_period_duration = mean(period_duration),
    .groups = "drop"
  )

For more information: https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#overdispersion

\[ NbBenthicDive \sim DayNight + MeanDiveDuration + MeanDayNightDuration + (1|ID) \]

Based on a table shape like the following:

ID	NbBenthic	DayNight	MeanDiveDuration	MeanDayNightDuration
A	10	Night	20	7
A	5	Day	21	17
B	11	Night	20	6
B	4	Day	22	18

Code

# model
mdl_50_extra <- glmer(
  nb_benthic_dive ~ day_night +
    scale(mean_dive_duration) +
    scale(mean_period_duration) +
    (1 | id),
  family = "poisson",
  data = data_plot %>% mutate(day_night = factor(day_night,
    # ordered = T,
    levels = c(
      "day",
      "night"
    )
  ))
)
# summary
broom.mixed::tidy(mdl_50_extra) %>%
  gt() %>%
  tab_header(title = "GLMM")

effect	group	term	estimate	std.error	statistic	p.value
GLMM
fixed	NA	(Intercept)	4.0982271	0.042556997	96.29972	0.000000e+00
fixed	NA	day_nightnight	-0.5263318	0.008752889	-60.13235	0.000000e+00
fixed	NA	scale(mean_dive_duration)	-0.1472508	0.011465763	-12.84265	9.457963e-38
fixed	NA	scale(mean_period_duration)	0.1262806	0.004555875	27.71820	4.213980e-169
ran_pars	id	sd__(Intercept)	0.8416426	NA	NA	NA

Code

qqnorm(residuals(mdl_50_extra))
qqline(residuals(mdl_50_extra))

Graph

Code

data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_100 == "Benthic") %>%
  mutate(
    hour = lubridate::hour(date),
    day = as.Date(date)
  ) %>%
  group_by(id, day, hour, day_night) %>%
  summarise(
    nb_benthic_dive = n(),
    mean_dive_duration = mean(Dduration),
    .groups = "drop"
  ) %>%
  group_by(id, day_night) %>%
  summarise(
    nb_average_benthic_dive = mean(nb_benthic_dive),
    mean_dive_duration = mean(mean_dive_duration),
    .groups = "drop"
  )


data_plot %>%
  ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +
  geom_boxplot(outlier.shape = NA) +
  labs(
    x = "Time of the day",
    y = "Average number of benthic dives per seal per hour"
  )

Figure 15: Boxplot of the number of (diff<100) benthic dives per hour and individual accross day-time and night-time

Model

Code

# determine what distribution would best describe our data
descdist(data_plot$nb_average_benthic_dive, discrete = FALSE)

summary statistics
------
min:  1   max:  5 
median:  2.565217 
mean:  2.612124 
estimated sd:  0.643816 
estimated skewness:  0.4024764 
estimated kurtosis:  3.504268

Figure 16: Result from `descdist` function to determine the best distribution to use for data modelling

Based on Figure 16, it seems that a gamma distribution would be the best to model this dataset.

\[ HourlyBenthicDiveRate \sim DayNight + MeanDiveDuration + (1|ID) \tag{3}\]

Based on a table shape like the following:

ID	HourlyBenthicDiveRate	DayNight	MeanDiveDuration
A	10	Night	20
A	5	Day	21
B	11	Night	20
B	4	Day	22

Code

# model
mdl_100 <- glmer(
  nb_average_benthic_dive ~ day_night +
    scale(mean_dive_duration) +
    (1 | id),
  family = Gamma,
  data = data_plot %>% mutate(day_night = factor(day_night,
    levels = c(
      "night",
      "day"
    )
  ))
)
# summary
broom.mixed::tidy(mdl_100) %>%
  gt() %>%
  tab_header(title = "GLMM")

effect	group	term	estimate	std.error	statistic	p.value
GLMM
fixed	NA	(Intercept)	0.41064106	0.003745655	109.631316	0.000000e+00
fixed	NA	day_nightday	-0.01556972	0.003443919	-4.520931	6.156812e-06
fixed	NA	scale(mean_dive_duration)	0.07981819	0.002733776	29.197048	2.113987e-187
ran_pars	id	sd__(Intercept)	0.03328018	NA	NA	NA
ran_pars	Residual	sd__Observation	0.13201902	NA	NA	NA

Go Further…

Code

data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_100 == "Benthic") %>%
  # I'm not so sure about this part..
  mutate(period_duration = if_else(
    day_night == "night",
    getNightLength(date, Lat, Long),
    24 - getNightLength(date, Lat, Long)
  )) %>%
  group_by(id, day_night) %>%
  summarise(
    nb_benthic_dive = n(),
    mean_dive_duration = mean(Dduration),
    mean_period_duration = mean(period_duration),
    .groups = "drop"
  )

For more information: https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#overdispersion

\[ NbBenthicDive \sim DayNight + MeanDiveDuration + MeanDayNightDuration + (1|ID) \]

Based on a table shape like the following:

ID	NbBenthic	DayNight	MeanDiveDuration	MeanDayNightDuration
A	10	Night	20	7
A	5	Day	21	17
B	11	Night	20	6
B	4	Day	22	18

Code

# model
mdl_100_extra <- glmer(
  nb_benthic_dive ~ day_night +
    scale(mean_dive_duration) +
    scale(mean_period_duration) +
    (1 | id),
  family = "poisson",
  data = data_plot %>% mutate(day_night = factor(day_night,
    # ordered = T,
    levels = c(
      "day",
      "night"
    )
  ))
)
# summary
broom.mixed::tidy(mdl_50_extra) %>%
  gt() %>%
  tab_header(title = "GLMM")

effect	group	term	estimate	std.error	statistic	p.value
GLMM
fixed	NA	(Intercept)	4.0982271	0.042556997	96.29972	0.000000e+00
fixed	NA	day_nightnight	-0.5263318	0.008752889	-60.13235	0.000000e+00
fixed	NA	scale(mean_dive_duration)	-0.1472508	0.011465763	-12.84265	9.457963e-38
fixed	NA	scale(mean_period_duration)	0.1262806	0.004555875	27.71820	4.213980e-169
ran_pars	id	sd__(Intercept)	0.8416426	NA	NA	NA

Code

qqnorm(residuals(mdl_100_extra))
qqline(residuals(mdl_100_extra))

Graph

Code

data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_150 == "Benthic") %>%
  mutate(
    hour = lubridate::hour(date),
    day = as.Date(date)
  ) %>%
  group_by(id, day, hour, day_night) %>%
  summarise(
    nb_benthic_dive = n(),
    mean_dive_duration = mean(Dduration),
    .groups = "drop"
  ) %>%
  group_by(id, day_night) %>%
  summarise(
    nb_average_benthic_dive = mean(nb_benthic_dive),
    mean_dive_duration = mean(mean_dive_duration),
    .groups = "drop"
  )


data_plot %>%
  ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +
  geom_boxplot(outlier.shape = NA) +
  labs(
    x = "Time of the day",
    y = "Average number of benthic dives per seal per hour"
  )

Figure 17: Boxplot of the number of (diff<150) benthic dives per hour and individual accross day-time and night-time

Model

Code

# determine what distribution would best describe our data
descdist(data_plot$nb_average_benthic_dive, discrete = FALSE)

summary statistics
------
min:  1   max:  5 
median:  2.565217 
mean:  2.610992 
estimated sd:  0.6399849 
estimated skewness:  0.4233921 
estimated kurtosis:  3.460216

Figure 18: Result from `descdist` function to determine the best distribution to use for data modelling

Based on Figure 18, it seems that a gamma distribution would be the best to model this dataset.

\[ HourlyBenthicDiveRate \sim DayNight + MeanDiveDuration + (1|ID) \tag{4}\]

Based on a table shape like the following:

ID	HourlyBenthicDiveRate	DayNight	MeanDiveDuration
A	10	Night	20
A	5	Day	21
B	11	Night	20
B	4	Day	22

Code

# model
mdl_150 <- glmer(
  nb_average_benthic_dive ~ day_night +
    scale(mean_dive_duration) +
    (1 | id),
  family = Gamma,
  data = data_plot %>% mutate(day_night = factor(day_night,
    levels = c(
      "night",
      "day"
    )
  ))
)
# summary
broom.mixed::tidy(mdl_150) %>%
  gt() %>%
  tab_header(title = "GLMM")

effect	group	term	estimate	std.error	statistic	p.value
GLMM
fixed	NA	(Intercept)	0.41135058	0.003639353	113.0285	0.000000e+00
fixed	NA	day_nightday	-0.01641792	0.003291418	-4.9881	6.097591e-07
fixed	NA	scale(mean_dive_duration)	0.08012749	0.002635001	30.4089	4.188961e-203
ran_pars	id	sd__(Intercept)	0.03289909	NA	NA	NA
ran_pars	Residual	sd__Observation	0.12728439	NA	NA	NA

Go Further…

Code

data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_150 == "Benthic") %>%
  # I'm not so sure about this part..
  mutate(period_duration = if_else(
    day_night == "night",
    getNightLength(date, Lat, Long),
    24 - getNightLength(date, Lat, Long)
  )) %>%
  group_by(id, day_night) %>%
  summarise(
    nb_benthic_dive = n(),
    mean_dive_duration = mean(Dduration),
    mean_period_duration = mean(period_duration),
    .groups = "drop"
  )

For more information: https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#overdispersion

\[ NbBenthicDive \sim DayNight + MeanDiveDuration + MeanDayNightDuration + (1|ID) \]

Based on a table shape like the following:

ID	NbBenthic	DayNight	MeanDiveDuration	MeanDayNightDuration
A	10	Night	20	7
A	5	Day	21	17
B	11	Night	20	6
B	4	Day	22	18

Code

# model
mdl_150_extra <- glmer(
  nb_benthic_dive ~ day_night +
    scale(mean_dive_duration) +
    scale(mean_period_duration) +
    (1 | id),
  family = "poisson",
  data = data_plot %>% mutate(day_night = factor(day_night,
    # ordered = T,
    levels = c(
      "day",
      "night"
    )
  ))
)
# summary
broom.mixed::tidy(mdl_150_extra) %>%
  gt() %>%
  tab_header(title = "GLMM")

effect	group	term	estimate	std.error	statistic	p.value
GLMM
fixed	NA	(Intercept)	4.1866975	0.042103917	99.43725	0.000000e+00
fixed	NA	day_nightnight	-0.5254813	0.008305595	-63.26835	0.000000e+00
fixed	NA	scale(mean_dive_duration)	-0.1149068	0.010908400	-10.53379	6.035410e-26
fixed	NA	scale(mean_period_duration)	0.1365982	0.004237414	32.23623	5.487327e-228
ran_pars	id	sd__(Intercept)	0.8334743	NA	NA	NA

Code

qqnorm(residuals(mdl_150_extra))
qqline(residuals(mdl_150_extra))

Graph

Code

data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_200 == "Benthic") %>%
  mutate(
    hour = lubridate::hour(date),
    day = as.Date(date)
  ) %>%
  group_by(id, day, hour, day_night) %>%
  summarise(
    nb_benthic_dive = n(),
    mean_dive_duration = mean(Dduration),
    .groups = "drop"
  ) %>%
  group_by(id, day_night) %>%
  summarise(
    nb_average_benthic_dive = mean(nb_benthic_dive),
    mean_dive_duration = mean(mean_dive_duration),
    .groups = "drop"
  )


data_plot %>%
  ggplot(aes(x = day_night, y = nb_average_benthic_dive)) +
  geom_boxplot(outlier.shape = NA) +
  labs(
    x = "Time of the day",
    y = "Average number of benthic dives per seal per hour"
  )

Figure 19: Boxplot of the number of (diff<200) benthic dives per hour and individual accross day-time and night-time

Model

Code

# determine what distribution would best describe our data
descdist(data_plot$nb_average_benthic_dive, discrete = FALSE)

summary statistics
------
min:  1   max:  5 
median:  2.541667 
mean:  2.600331 
estimated sd:  0.6434467 
estimated skewness:  0.3851526 
estimated kurtosis:  3.429033

Figure 20: Result from `descdist` function to determine the best distribution to use for data modelling

Based on Figure 20, it seems that a gamma distribution would be the best to model this dataset.

\[ HourlyBenthicDiveRate \sim DayNight + MeanDiveDuration + (1|ID) \tag{5}\]

Based on a table shape like the following:

ID	HourlyBenthicDiveRate	DayNight	MeanDiveDuration
A	10	Night	20
A	5	Day	21
B	11	Night	20
B	4	Day	22

Code

# model
mdl_200 <- glmer(
  nb_average_benthic_dive ~ day_night +
    scale(mean_dive_duration) +
    (1 | id),
  family = Gamma,
  data = data_plot %>% mutate(day_night = factor(day_night,
    levels = c(
      "night",
      "day"
    )
  ))
)
# summary
broom.mixed::tidy(mdl_200) %>%
  gt() %>%
  tab_header(title = "GLMM")

effect	group	term	estimate	std.error	statistic	p.value
GLMM
fixed	NA	(Intercept)	0.41345050	0.003696624	111.845416	0.000000e+00
fixed	NA	day_nightday	-0.01702406	0.003423988	-4.971998	6.626632e-07
fixed	NA	scale(mean_dive_duration)	0.08153323	0.002723215	29.940067	5.926044e-197
ran_pars	id	sd__(Intercept)	0.03261119	NA	NA	NA
ran_pars	Residual	sd__Observation	0.13047509	NA	NA	NA

Go Further…

Code

data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_200 == "Benthic") %>%
  # I'm not so sure about this part..
  mutate(period_duration = if_else(
    day_night == "night",
    getNightLength(date, Lat, Long),
    24 - getNightLength(date, Lat, Long)
  )) %>%
  group_by(id, day_night) %>%
  summarise(
    nb_benthic_dive = n(),
    mean_dive_duration = mean(Dduration),
    mean_period_duration = mean(period_duration),
    .groups = "drop"
  )

For more information: https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#overdispersion

\[ NbBenthicDive \sim DayNight + MeanDiveDuration + MeanDayNightDuration + (1|ID) \]

Based on a table shape like the following:

ID	NbBenthic	DayNight	MeanDiveDuration	MeanDayNightDuration
A	10	Night	20	7
A	5	Day	21	17
B	11	Night	20	6
B	4	Day	22	18

Code

# model
mdl_200_extra <- glmer(
  nb_benthic_dive ~ day_night +
    scale(mean_dive_duration) +
    scale(mean_period_duration) +
    (1 | id),
  family = "poisson",
  data = data_plot %>% mutate(day_night = factor(day_night,
    # ordered = T,
    levels = c(
      "day",
      "night"
    )
  ))
)
# summary
broom.mixed::tidy(mdl_200_extra) %>%
  gt() %>%
  tab_header(title = "GLMM")

effect	group	term	estimate	std.error	statistic	p.value
GLMM
fixed	NA	(Intercept)	4.2052051	0.042347958	99.30125	0.000000e+00
fixed	NA	day_nightnight	-0.5284598	0.008190191	-64.52351	0.000000e+00
fixed	NA	scale(mean_dive_duration)	-0.1162468	0.010784403	-10.77916	4.318092e-27
fixed	NA	scale(mean_period_duration)	0.1405915	0.004144919	33.91899	3.496541e-252
ran_pars	id	sd__(Intercept)	0.8386222	NA	NA	NA

Code

qqnorm(residuals(mdl_200_extra))
qqline(residuals(mdl_200_extra))

7 Boxplot of the number of benthic dive (day vs. night)

Code

data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName == "Benthic") %>%
  group_by(id, day_night) %>%
  summarise(nb_benthic_dive = n(), .groups = "drop")

data_plot %>%
  ggplot(aes(x = day_night, y = nb_benthic_dive)) +
  geom_boxplot(outlier.shape = NA) +
  coord_cartesian(ylim = c(0, 1000)) +
  labs(
    x = "Time of the day",
    y = "Number of benthic dives per individual"
  )

Figure 21: Boxplot of the average number of (all) benthic dives per hour and individual accross day-time and night-time

Code

broom::tidy(glm(
  nb_benthic_dive ~ day_night,
  family = "quasipoisson",
  data = data_plot
)) %>%
  gt() %>%
  fmt_number(
    columns = 2:5,
    decimals = 2,
    suffixing = TRUE
  ) %>%
  tab_header(title = "nb_benthic_dive ~ day_night")

term	estimate	std.error	statistic	p.value
nb_benthic_dive ~ day_night
(Intercept)	5.96	0.05	117.16	0.00
day_nightnight	−0.84	0.09	−9.07	0.00

Code

data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_50 == "Benthic") %>%
  group_by(id, day_night) %>%
  summarise(nb_benthic_dive = n(), .groups = "drop")

data_plot %>%
  ggplot(aes(x = day_night, y = nb_benthic_dive)) +
  geom_boxplot(outlier.shape = NA) +
  coord_cartesian(ylim = c(0, 250)) +
  labs(
    x = "Time of the day",
    y = "Number of benthic dives per individual"
  )

Figure 22: Boxplot of the average number of (diff<50) benthic dives per hour and individual accross day-time and night-time

Code

broom::tidy(glm(
  nb_benthic_dive ~ day_night,
  family = "quasipoisson",
  data = data_plot
)) %>%
  gt() %>%
  fmt_number(
    columns = 2:5,
    decimals = 2,
    suffixing = TRUE
  ) %>%
  tab_header(title = "nb_benthic_dive ~ day_night")

term	estimate	std.error	statistic	p.value
nb_benthic_dive ~ day_night
(Intercept)	4.54	0.07	65.08	0.00
day_nightnight	−0.54	0.12	−4.69	0.00

Code

data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_100 == "Benthic") %>%
  group_by(id, day_night) %>%
  summarise(nb_benthic_dive = n(), .groups = "drop")

data_plot %>%
  ggplot(aes(x = day_night, y = nb_benthic_dive)) +
  geom_boxplot(outlier.shape = NA) +
  coord_cartesian(ylim = c(0, 250)) +
  labs(
    x = "Time of the day",
    y = "Number of benthic dives per individual"
  )

Figure 23: Boxplot of the average number of (diff<100) benthic dives per hour and individual accross day-time and night-time

Code

broom::tidy(glm(
  nb_benthic_dive ~ day_night,
  family = "quasipoisson",
  data = data_plot
)) %>%
  gt() %>%
  fmt_number(
    columns = 2:5,
    decimals = 2,
    suffixing = TRUE
  ) %>%
  tab_header(title = "nb_benthic_dive ~ day_night")

term	estimate	std.error	statistic	p.value
nb_benthic_dive ~ day_night
(Intercept)	4.61	0.07	65.55	0.00
day_nightnight	−0.54	0.12	−4.65	0.00

Code

data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_150 == "Benthic") %>%
  group_by(id, day_night) %>%
  summarise(nb_benthic_dive = n(), .groups = "drop")

data_plot %>%
  ggplot(aes(x = day_night, y = nb_benthic_dive)) +
  geom_boxplot(outlier.shape = NA) +
  coord_cartesian(ylim = c(0, 250)) +
  labs(
    x = "Time of the day",
    y = "Number of benthic dives per individual"
  )

Figure 24: Boxplot of the average number of (diff<150) benthic dives per hour and individual accross day-time and night-time

Code

broom::tidy(glm(
  nb_benthic_dive ~ day_night,
  family = "quasipoisson",
  data = data_plot
)) %>%
  gt() %>%
  fmt_number(
    columns = 2:5,
    decimals = 2,
    suffixing = TRUE
  ) %>%
  tab_header(title = "nb_benthic_dive ~ day_night")

term	estimate	std.error	statistic	p.value
nb_benthic_dive ~ day_night
(Intercept)	4.64	0.07	65.52	0.00
day_nightnight	−0.54	0.12	−4.63	0.00

Code

data_plot <- data_dive %>%
  # get rid of data without location information
  filter(!is.na(Lat)) %>%
  # keep all benthic dives
  filter(DiveTypeName_200 == "Benthic") %>%
  group_by(id, day_night) %>%
  summarise(nb_benthic_dive = n(), .groups = "drop")

data_plot %>%
  ggplot(aes(x = day_night, y = nb_benthic_dive)) +
  geom_boxplot(outlier.shape = NA) +
  coord_cartesian(ylim = c(0, 250)) +
  labs(
    x = "Time of the day",
    y = "Number of benthic dives per individual"
  )

Figure 25: Boxplot of the average number of (diff<200) benthic dives per hour and individual accross day-time and night-time

Code

broom::tidy(glm(
  nb_benthic_dive ~ day_night,
  family = "quasipoisson",
  data = data_plot
)) %>%
  gt() %>%
  fmt_number(
    columns = 2:5,
    decimals = 2,
    suffixing = TRUE
  ) %>%
  tab_header(title = "nb_benthic_dive ~ day_night")

term	estimate	std.error	statistic	p.value
nb_benthic_dive ~ day_night
(Intercept)	4.67	0.07	65.47	0.00
day_nightnight	−0.55	0.12	−4.66	0.00

8 Extra

Let’s look at the distribution of the dive duration to decide a threshold.

Code

data_dive %>%
  ggplot(aes(x = Dduration)) +
  geom_histogram(
    fill = "grey80",
    color = "black"
  ) +
  geom_vline(xintercept = 23 * 60, linetype = "dashed", col = "lightgreen") +
  geom_vline(xintercept = 4000, linetype = "dashed", col = "salmon") +
  labs(y = "# of dives", x = "Dive duration (s)") +
  theme_minimal() +
  theme(axis.line = element_line(arrow = arrow(type = "closed", length = unit(0.1, "inches"))))

Ok, so let’s pick 2500 sec and look on a map where this dives are.

Code

# define palette
data_inter <- data_dive %>%
  filter(!is.na(Lat) &
    Dduration > 4000) %>%
  arrange(id, date) %>%
  mutate(
    date_post = date + Dduration,
    diff_time = as.numeric(
      c(difftime(
        tail(date, -1),
        head(date_post, -1)
      ), 0),
      units = "mins"
    )
  )

colPal <- scales::col_numeric(
  palette = "RdYlGn",
  domain = c(0, max(data_inter$diff_time))
)
trip +
  # new scale
  new_scale_fill() +
  # kernel
  ggspatial::geom_spatial_point(
    data = data_inter,
    aes(x = Long, y = Lat),
    crs = 4326,
    col = colPal(data_inter$diff_time),
    size = 1
  ) +
  # no display alpha
  guides() +
  # title
  labs(
    title = "Dives longer than 4000 seconds",
    subtitle = paste0(nrow(data_inter), " dives")
  )

1 Import library

2 Setting up custom function

2.1 windrose

2.2 getNightLength

3 Import data

4 Maps

5 Windrose

6 Boxplot of the benthic dive rate (day vs. night)

Graph

Model

Go Further…

Graph

Model

Go Further…

Graph

Model

Go Further…

Graph

Model

Go Further…

Graph

Model

Go Further…

7 Boxplot of the number of benthic dive (day vs. night)

8 Extra

2.1 `windrose`

2.2 `getNightLength`